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ABSTRACT 



Context. The presence of vortices in accretion discs has been debated for more than a decade. Baroclinic instabilities might be a way 
to generate these vortices in the presence of a radial entropy gradient. However, the nature of these instabilities is still unclear and 3D 
parametric instabilities can lead to the rapid destruction of these vortices. 

Aims. We present new results exhibiting a subcritical baroclinic instability (SBI) in local shearing box models. We describe the 2D 
and 3D behaviour of this instability using numerical simulations and we present a simple analytical model describing the underlying 
physical process. 

Methods. We investigate the SBI in local shearing boxes, using either the incompressible Boussinesq approximation or a fully com- 
pressible model. We explore the parameter space varying several local dimensionless parameters and we isolate the regime relevant 
for the SBI. 3D shearing boxes are also investigated using high resolution spectral methods to resolve both the SBI and 3D parametric 
instabilities. 

Results. A subcritical baroclinic instability is observed in flows stable for the Solberg-Hoi'land criterion using local simulations. This 
instability is found to be a nonlinear (or subcritical) instability, which cannot be described by ordinary linear approaches. It requires a 
radial entropy gradient weakly unstable for the Schwartzchild criterion and a strong thermal diffusivity (or equivalently a short cooling 
time). In compressible simulations, the instability produces density waves which transport angular momentum outward with typically 
a < 3 x 10~ 3 , the exact value depending on the background temperature profile. Finally, the instability survives in 3D, vortex cores 
becoming turbulent due to parametric instabilities. 

Conclusions. The subcritical baroclinic instability is a robust phenomenon, which can be captured using local simulations. The in- 
stability survives in 3D thanks to a balance between the 2D SBI and 3D parametric instabilities. Finally, this instability can lead to a 
weak outward transport of angular momentum, due to the generation of density waves by the vortices. 
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1. Introduction 

The existence of long-lived vort i ces in accretion discs was 
first proposed by Ivon Weizsackerl (Il944l) i n a model of planet 
form ation. This idea was reintroduced by IBarge & Sommerial 
(Il995l) to accelerate the formation of planetesimals through 
a dust trapping process. Moreover, large scale vortices can 
lead to a significant transport of angul ar momentum thanks t o 
the production of density waves dJohnson & Gammiell20 05b). 
Many physical processes have been introduced in the lit- 
erature to justify the existe nce of these vortic es such as 
Rossby wave in stabilities (lLovelace et alJ Il999l) . planetary 
gap instabilitie s (Ide V al-Borr o etaD l2006l I2007I) . 3D circu- 
lation models (iBarranco & Marcusl l2005l) . MHD turbulence 
(Fromang & Nelson 2005) and the global baroclinic instability 
( Klahr & Bodenheime ri2003l) . 

Baroclinic instabilities in accretion discs has regained in- 
terest in the past few years. A first version of these instabil- 
ities (the global baroclinic insta bility or GBI) was mentioned 
by iKlahr & Bodenheimerl (I2003I) using a purely numerical ap- 
proach and considering a disc with a radial entropy gradient. 
However, the presence of this instability was affected by chang- 
ing the numerical scheme used in the simulations, casting strong 
doubts on this instability as a real physical proc esses. The li near 
properties of the GBI were th en investigated bv lKlahrl(l2004 and 
Johnson & Gammie ( 2005a). However, only transient growths 



were found in this context. iKlahrl ([2004) speculated that these 
transient growt hs could lead to a nonlinear in stability explain- 
ing the result of Klahr & Bodenheimer] (I2003I) . but it is still un- 
clear whether tr ansient growths are releva nt for nonlinear insta- 
bilities (see e.g. iLesur & Lon garetti 2005, and reference therein 
for a complete discussion of this point). T he nonlinear problem 
was then studied in local shearing boxes by I Johnson & Gam mie 
(2006) using fully compressible numerical methods. These au- 
thors found no instability in the Keplerian disc regime and con- 
cluded that the GBI was "either global or nonexistent". 

Neverthele s s, baro clinic processes were studied again by 
iPetersen et al.l d2007a) using anelastic global simulations and 
includ ing new physical processes. As in Klahr & Bodenheimer 
( 2003), a weak radial entropy gradient was imposed in the sim- 
ulation. However, a cooling function was also included in their 
model, in order to force the system to rel ax to the imposed radia l 
temperature profile. In these simulations. IPetersen et al.l (l2007ah 
observed spontaneous formation of vortices with radial entropy 
gradients compatibl e with accretion disc thermodynamics. In 
a subsequent paper (Pet ersen et al.ll2007bb . these vortices were 
found to survive for several hu ndreds orbits. According to the 
authors, their disagreement with lJohnson & Gam mie ( 2006|) was 
due to a larger Reynolds number a nd the use of spectral method s, 
a possibility already mentioned by Johnson & Gammie (I2006I) . 

In this paper, we revisit baroclinic instabilities in a local 
setup (namely the shearing box) both in the incompressible- 
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Boussinesq approximation and in fully compressible simula- 
tions. The aim of this paper is to clarify several of the points 
which have been discussed in the literature for the past 6 years 
about the existence, the nature and the properties of baroclinic 
instabilities. We show that a subcritical baroclinic instability (or 
shortly SBI) does exist in shearing boxes. This instability is non- 
linear (or subcritical) and stro ngly linked to t herma l diffusion, 
a point already mentioned by iPetersen et al.l (l2007bl) . We also 
present results concerning the behaviour of the SBI in 3 dimen- 
sions, as vortice s are known to be unstable because of paramet- 
ric instabilities dLesur & Papaloizoull2009l) . We begin in §[2] by 
presenting the equations, introducing the dimensionless num- 
bers and describing the numerical methods used in the rest of 
the paper. §|3]is dedicated to 2D incompressible-Boussinesq sim- 
ulations and a qualitative understanding of the instability. We 
present in §|4] our results in compressible shearing boxes and in 
§|5]the 3D behaviour of the instability. Conclusions and implica- 
tions of our findings are discussed in §0 

2. Local model 

2.1. Equations 

In the following, we will assume a local model for the ac- 
cretion disc, followi ng the shearing shee t approximation . The 
reader may consul t Hawlev et al.l (Il995l) . iBalbusI (I2003I) and 
Reg ev & Umurhanl ([2008) for an extensive discussion of the 
properties and limitations of this model. As a simplification, 
we will assume the flow is in compressible, consistently with 
the small shearing box model ( Re gev & Umurhanl l2008h . The 
shearing-box equations are found by considering a Cartesian 
box centred at r = Rq, rotating with the disc at angular velocity 
D = £l(Ro). We finally introduce in this box a radial stratificatio n 
using the Boussinesq approximation (Spiegel & Veronisl 1960b . 
Defining r - Rq — > x and Ro(p — > y as in Haw lev et al.l ( 19951) . 
one eventually obtains the following set of equations: 



d t u + u- Vu = - VII - 2Q x u 

2£lSxe x - AN 2 6e x + vAu, 
d t 6 + u-V6 = u x /A + fiAO 
V • u = 0, 



(1) 

(2) 
(3) 



where u = u x e x +u y e y +u z e z is the fluid velocity, 6 the potential 
temperature deviation, v the kinematic viscosity and [i the ther- 
mal diffusivity. In these equations, we have defined the mean 
shear S = -r<9 r Q, which is set to S - (3/2)11 for a Keplerian 
disc. One can check easily that the velocity field U = -S xe y 
is a steady solution of these equations. In the following we will 
consider the evolutions of the perturbations v (not necessarily 
small) of this profile defined by v = u - U. 

The generalised pressure II is calculated solving a Poisson 
equation with the incompressibility condition ([3]). For homo- 
geneity and consistency with the traditional Boussinesq ap- 
proach, we have introduced a stratification length A. Note how- 
ever that A disappears from the dynamical properties of these 
equations as one can renormalise the variables defining ff = AO. 
The stratification itself is controlled by the Brunt- Vaisala fre- 
quency N, defined for a perfect gas by 



AT 



I dP d /P\ 



P 



(4) 



where P and p are assumed to be the background equilibrium 
profile and y is the adiabatic index. With these notations, one 
can recover the ordinary thermodynamical variables as 6 = dp/p, 



where dp is the density perturbation and p is the background den- 
sity. The stratification length is then defined by A = -OrP/pN 2 . 

2.2. Dimensionless numbers 

The system described above involves several physical processes. 
To clarify the regime in which we are working and to make eas- 
ier comparisons with previous work, we define the following di- 
mensionless numbers: 

- The Richardson number Ri = N 2 /S 2 compares the shear 
timescale to the buoyan cy timescale. This definition i s equiv- 
alent to the definition of Jo hnson & Gammiel (l2005al) . 

- The Peclet number Pe = SL 2 /fi. 

- The Reynolds number Re = SL 2 /v. 

L is a typical scale of the system, chosen to be the radial box size 
(L x ) in our notations. 

The linear stability properties of this flow are quite well un- 
derstood. The flow is linearly stable for axisymmetric perturba- 
tions when the Solberg-Hoiland criterion is satisfied. This crite- 
rion may be written locally as: 

2£l(2Q-S) + N 2 > 



Stability, 



(5) 



or equivalently in a Keplerian disc, Ri > -4/9. Another linear 
stability criterion, the Schwarzschild criterion, is often used for 
convectively stable flows without rotation nor shear. This crite- 
rion reads 



N 2 > Stability, 



(6) 



or equivalently Ri > 0. As we will see in the following, this 
criterion is the relevant one for the SBI. 

When viscosity and thermal diffusivity are included, 
the Solberg-Hoiland criterion is modified, potentially lead- 
ing to the Goldrei c h-Sch u bert-Fri c ke (G SF) instability 
dGoldreich & Schubertl 1 19671: iFrickd 119681) . The stability 
criterion is in that case 

yu2Q(2Q - S ) + vN 2 > Stability. (7) 

This criterion is satisfied when both © is verified and v/fi < 1, 
which corresponds to the regime studied in this paper. 

2.3. Numerical methods 

We have used two different codes to study this instability. 
When using the Boussineq model (O and §0) , we have used 
SNOOPY, a 3D incompressible spectral code. This code uses 
an implicit scheme for thermal and viscous diffusion, allowing 
us to study simulations with small Pe without any strong con- 
strain on the CFL condition. The spectral algorithm of SNOOPY 
has now been used in s everal hydrodynamic and magn etohy- 
drodynamic studies (e.g. iLesur & Longarettill2005l 12007) and is 
freely available on the author's website . For compressible sim - 
ulations (©, we have used NIRVANA dZiegler& Yorkdl 1997b . 
NIRVANA has been used frequently in the past to study var- 
ious problems involving MHD turbulence in the she aring box 
dFromang & Papaloizoull2006t iPapaloizou et al.ll2004l) . 

In t he following, we u se the shearing sheet boundary condi- 
tions (Hawlev et al. 1995) in the radial direction. This is made 
possible by the use of the Boussinesq approximation in which 
only the gradients of the background profile appear explicitly, 
through A and N. Therefore, when A and N are constant through 
the box, one can assume that the thermodynamic fluctuation is 
shear-periodic, consistently with the shearing- sheet approxima- 
tion. 
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Fig. 1. Volume averaged enstrophy for several initial amplitude 
perturbation (Ap) in arbitrary units. A subcritical instability is 
observed for finite amplitude perturbations between Ap = 0.2 
and Ap = 0.4. 

3. 2D subcritical baroclinic instability in 
incompressible flows 

To start our investigation, we consider the simplest case of a 
2D (x, y) problem in an infinitely thin disc. This s etup is the lo- 
cal eq uivalent of the 2D global anelastic setup of iPetersen et alJ 
(l2007al) . and one expects to find similar properties in both cases 
if the instability is local. In two dimensions, it is easier to con- 
sider the vorticity equation instead of (Q]). Defining the vertical 
vorticity of the pertubations by uj = d x v y - d y v x , we have: 

d t D + vV(i)-S xdyoj = AN 2 d y 6 + vAoj (8) 

In this formulation, the only source of enstrophy (oj 2 /2) = 
J dxdyaj 2 /2 is the baroclinic term AN 2 d y 6 which will be shown 
to play an important role for the instability. 

The simulations presented in this section are computed in a 
square domain L x = L y = 1 with a resolution of512x512 grid 
points using our spectral code. When not explicitly mentioned, 
the time unit is the shear timescale S~ l . We choose the fiducial 
parameters to be Re = 4 x 10 5 ; Pe = 4 x 10 3 and Ri = -0.01, 
in order to have a flow linearly stable for the Solberg-Hoi'land 
criterion ([5]). These p arameters are close to the one used by 
IPetersen et al.l d2007a) and are compatible with disc thermody- 
namics. 

3.1. Influence of initial perturbations 

In the first numerical experiment, we choose to test the ef- 
fect of the initial conditions, keeping constant the dimension- 
less parameters. In our initial conditions, we excite randomly 
the largest wavelengths of the vorticity field with an amplitude 
A p. This initial condition is slightly different from the one used 
by IPetersen et al.l (l2007allbh who introduced perturbations in the 
temperature field. In each experiment we modify the amplitude 
of the initial perturbation Ap, and follow the time evolution of 
the total enstrophy (Fig.[T]) 

The numerical results indicates clearly the presence of a 
nonlinear or subcritical transition in the flow. Indeed, we find 
the "instability" for lar ge enough initial perturbatio ns. This also 
confirms the result of Johnson & Gammie (I2006I) : there is no 



linear instability in the presence of a weak radial stratification . 
Moreover, one of the reasons that iJohnson & Gammiel (I2006I) 
did not observe any transition could be because of the weak 
initial perturbations they used (between 10" 12 and 10" 4 ). Using 
similar initial conditions, we did not observe any transition ei- 
ther. The existence of the instability in our simulations was 
checked by doubling the resolution (1024 x 1024), keeping con- 
stant all the dimensionless numbers. We found no significant dif- 
ference between high and low resolution results, showing that 
our simulations were converged for this problem. 

When the system undergoes a subcritical transition, it de- 
velops long-lived self-sustained vortices, as shown for Ap = 1 
in Fig. [2] (top). For such a large Reynolds n umber, it is known 
that v ortices survive for long times (see e.g. lUmur han & Regevl 
12004 . even without any baroclinic effect. To check that the ob- 
served vortices were really due to the baroclinic term, we have 
carried out the exact same simulation but without stratification 
(Fig.[2]bottom). This simulation also shows the formation of vor- 
tices but these are lately dissipated on a few hundred shear times. 
At t = 500 S~ l the difference between the simulation with and 
without baroclinicity becomes obvious. 

The vortices observed in these incompressible simulations 
lead to a very weak and strongly oscillating turbulent transport of 
angular momentum. One finds typically an inward transport with 
a = (v x v v )/S L 2 ^ -3 x 10~ 5 . This is consistent with the results 
published by IPetersen et al.l (l2007bl) in global simulations. 

In several other numerical experiments (not shown here), we 
have noticed that the amplitude threshold required to get the in- 
stability depends on Re and Pe, a larger Reynolds number being 
associated with a weaker initial perturbati on. This dependency 
was pointed out by IPetersen et al.l (l2007al) and it indicates that 
the amplitude threshold in a realistic disc could be very small 
(i.e. much smaller than the sound speed). Note however that this 
threshold might also be scale dependent, a problem which has 
not been addressed here. 

3.2. Influence of the Richardson number 

To understand how the instability depends on the amplitude 
of the baroclinic term, we have carried out several runs with 
Re = 4 x 10 5 and Pe = 4 x 10 3 , varying Ri from -0.02 to 
-0.16 and from 0.02 to 0.016. We compare the resulting shell- 
integrated enstrophy spectra (<^/2) on Fig. [4] When Ri < and 
\Ri\ becomes larger, the instability tends to amplify smaller and 
smaller scales. In particular for Ri = -0.02, the dominant scale 
is clearly the box scale, as already observed for the fiducial run. 
This trend is also observed looking directly at vorticity snapshots 
(Fig. [3]). The sign of Ri is also of importance for the onset of the 
SBI. To demonstrate this effect, we show the time history of the 
total enstrophy for positive and negative Ri on figure In the 
cases of positive Ri, the total enstrophy decays until the flow be- 
comes axisymmetric. Perturbations are then damped very slowly 
on a viscous timescale. We conclude from these results that a 
necessary condition for the SBI to appear is a flow which do not 
satisfy the Schwarzschild criterion ©. 

3.3. Influence of thermal diffusion 

The importance of th ermal cooling and ther mal diffusion was 
already pointed out by IPetersen et al.l (l2007bl) . To check this de- 
pendancy in our local model, we have considered several simula- 
tions with Ri = -0.01, varying Pe from Pe = 20 to Pe = 16000. 
The resulting enstrophy evolution is presented for several of 
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Fig. 2. Evolution of the vorticity in the fiducial case. Top row has a baroclinic term with Ri = -0.01 . Bottom row has no baroclinicity. 
We show t = 10 (left), t = 100 (middle), t = 500 (right). 
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Fig. 3. Vorticity map for Ri = -0.02 (left) and Ri = -0.16 (right) at t = 500. As already shown by the spectra, small scales are 
dominant for larger \Ri\. 



these runs on Fig. [6j Looking at the snapshots of the vorticity 
field for these simulations, we find the SB I approximatively for 
50 < Pe < 8000. Assuming a typical vortex size / ~ 0.25 (see 
e.g. Fig 12), these Peclet numbers correspond to thermal diffusion 
times 3 S~ l < r^s < 500 S~ l over the vortex size, with an opti- 
mum found f orrdiff - IPS" 1 {Pe = 250). Note that the cooling 
time used by iPetersen et al.l d2007bl) lies typically in this range 
of values. 



3.4. Instability mechanism 

Since the flow is subcritically unstable, no linear anal- 
ysis can capture this in stability entirely. As shown by 



plification going asymptotically lik43 \Ri\t l ~ 4Rl for \Ri\ «: 1 
when ji = v = 0, the waves being ultimately decaying when 
ji > or y > 0. However, this tells us little to nothing about 
the nonlinear behaviour of such a flow. Indeed, it is known that 
barotropic Keplerian fl ows undergo arbitrarily l arge transient 
amplifications (see e.g. IChagelishvili et al.l 120031), but subcrit- 
ical transitions are yet to be found in that ca se (lHawley et al.l 
119991: iLesur & Longarettill2005l: iJi et al.ll2QQ6b . In the SBI case, 
the subcritical transition happens only for negative Ri as shown 
above, in fl agrant contradiction with th e linear amplification de- 
scribed by John son & Gammiel (l2005ah . This is a strong indica- 
tion that linear theory is of little use for describing this instabil- 
ity. 

To isolate a mechanism for a subcritical instability, one 
should start from a non-linear structure observed in the simu- 



Joh nson & Gammiel (l2005al) . an ensemble of shearing waves in 
a baroclinic flow is subject to a transient growth with an am- 



1 Note that the ampli fication occurs independently of the sign of Ri, 
as already mentioned by Johns on & Gam mie ( 2005a|). 
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Fig. 4. Enstrophy spectrum as a function of the Richardson num- Fig. 6. Time history of the total enstrophy for several Peclet 
ber RL As \Ri\ becomes larger, the instability moves to smaller number (Pe). The instability is observed for 50 < Pe < 8000. 
scales. 
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Fig. 5. Time history of the total enstrophy for positive and neg- 
ative Richardson number. The instability is observed only for 
negative Ri. 

lations. For the SBI, these structures are self- sustained vortices. 
In order to understand how baroclinic effects can feedback on 
the vortex structure, we initi alised our fi ducial simulation with a 
Kida vortex of aspect ratio 4 (Kid a!l98ll) . Although this vortex is 
known to be an exact nonlinear solution of the inviscid equation 
of motions, it is slowly modified by the explicit viscosity and the 
baroclinic term, leading to a slow growth of the vortex. We show 
on Fig. |7] the resulting vortex structure (left) and the associated 
baroclinic term AN 2 d y (right) at t = 100 S" 1 . Note that these 
structures are quasi steady and evolve on timescales much longer 
than the shear timescale. It is clear from these snapshots that the 
baroclinic feedback tends to amplify the vorticity located inside 
the vortex, leading to the growth of the vortex itself. 

To understand the origin of the baroclinic feedback, let's as- 
sume the physical quantities are evolving slowly in time, as is 
observed in the simulations, so we may write 



oj = co(et) = oj(t), 
6 = 6(et) = 6(t), 



(9) 
(10) 



where 6 is a small parameter and it is assumed that \6 T \ is of the 
same order as |v • V0|, and similarly for to. We therefore have to 
solve: 



ed T oj + (v • V - Sxd y )oj = AN d y 6 + vAqj, 
ed T 6 + (v • V - Sxd y )0 = v x /A + jiAO. 



(11) 
(12) 



Since the Richardson number is assumed to be small, the baro- 
clinic feedback of (Tf2t has to be small. As only this term can 
lead to an instability, we can assume it scales like e. The role of 
the viscosity is to prevent the instability from happening, damp- 
ing vorticity fluctuations. Assuming we are in a regime in which 
the instability appears, the viscosity has to be of the order of the 
baroclinic term, scaling like e. At zeroth order in 6, we are left 
with: 



(v-V-Sxd y )cD = 0, 
(v • V - Sxd y )0 = v x /A + jiAO 



(13) 
(14) 



This system of equations describes a quite simple physical sys- 
tem: the vorticity field is a steady solution of the inviscid vor- 
ticity equation with a constant shear and without any baro- 
clinic feedback, whereas the potential temperature results from 
advection-diffusion of the background entropy profile due to the 
flow structure. A family of steady solutions of the inviscid vor- 
ticity equation with a consta nt shear are known to be steady 
vortices, like the iKidal ([1981) vortex. More generally, it can be 
shown that any cd of the form: 



co = F(il/) 



(15) 



where if/ is the stream function defined by u = e z X is 
solution of (fT3l) . 

In the following, we assume (TT3l) is satisfied by co and we 
look for a solution to the entropy equation (fT4l) . As a further sim- 
plification, we neglect the advection of the perturbed entropy 6 
by v, assuming this effect is small compared to the advection by 
the mean shear, although this does not affect the argument lead- 
ing to the possibility of instability (see below). Using a Fourier 
decomposition 6{x) = J dk 6(k) exp(//c • cc), one gets 



Sk y —+Lik 1 6 = v x /A. 
dk x 



(16) 
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Fig. 7. Vortex structure obtained starting from a Kida vortex, (left) vorticity map. (right) baroclinic term (AN 2 d y 6) map. 



A solution to this equation can easily be found using standard 
techniques for solving first order ordinary differential equations. 
Assuming |v x | decays to zero when \k\ —> oo, one may write the 
result in the form 



0(k x ,k y ) = 



1 



ex 



dk' x V X (K,ky)G(k X ,k' x ,ky), 



A\Sk y \ 
where we have defined 
G(k x ,k' x ,k y ) = 9i((k x -k' x )Sk y /\Sk y \y, 



(17) 



exp 



(18) 



H being the Heaviside function. From this expression, it is for- 
mally possible to derive the time evolution of the vorticity and 
look for an instability using the first order term in (IT2b . However, 
this requires an explicit expression for the vorticity, which is a 
priori unknown. Instead, we have choosen to compute the evolu- 
tion of the total enstrophy of the system (oj 2 /2) = J dxdy oj 2 /2: 



d t (a) 2 /2) = KN 2 {cod y 6) - v<|VM > 



(19) 



Evidently, an instability can occur only if the first term on the 
right hand side is positive. Using (fTTlh is is possible to derive an 
analytical expression for this term: 



AN 2 (a)d y 6) = -4n 2 % 
\k y \ 



dk x dk y dk x 



\S\\k 



'|2 



cb*(k)G(k x ,k x ,k y )a)(k') 



(20) 



where k f - (k x , k y ). In this expression, the existence of an in- 
stability is controlled by the Green's function G(k x , k x ,k y ) and 
the shape of a>(k). To our knowledge, it is not possible to de- 
rive any general criterion for the instability, as such a criterion 
would depend on the background vortex considered. It is how- 
ever possible to derive a criterion in the limit of a large thermal 
diffusivity \i. In this limit, G is strongly peaked around k' x - k x , 
and one can assume in first approximation that G is proportional 
to a Dirac distribution. More formally, one can approximate in 
the limit of large thermal diffusivity: 



G(k x ,k' x ,k y ) - <n((k x - k' x )Sk y /\sk y \yx V 



Jk k2(K ~ kx) 



(21) 



when k x k' x . It is then possible to derive an expression for the 
baroclinic feedback as a series expansion in the small parameter 



r C i k 

LN 2 {cod y 6) * -A7i 2 %yN 2 j dk \^4 

Ski 



^k 4 



&(k)\ 2 - 

*{k)d kx [cb(k)lk 2 ] + ...} 



(22) 



In this approximation, we find two competing effects. The first 
term on the right hand side appears when the thermal diffusion 
completely dominates over the shear in ([T6l) . It has a positive 
feedback on the total enstrophy when N 2 < and can be seen as 
a source term for the SBI. The second term involves a competi- 
tion between shear and diffusion. It can be seen in first approx- 
imation as a phase mixing term, and the resulting can be either 
positive or negative, depending on the vorticity background one 
considers. Physically, this term shears out the entropy structure 
created by the vortex and if it is too strong, it kills the positive 
feedback of the first term. One should note however that this 
term will not reverse the sign of the baroclinic feedback but will 
just weaken it by randomising the phase coherence between a> 
and 0. 

To understand the "growth" described by (l22b , on can de- 
rive an order of magnitude expression for the growth rate y = 
d In (oj 2 )/dt combining ([191) with (l22l) : 



r 



(-N 2 )o- 2 



(S a 2 //i)-—. 



(23) 



In this expression, cr is the typical vortex size, with the assump- 
tion cr ~ l/k. We have included a phase mixing term with the 
function as an extension of (l22t . This function depends on 
the background vortex solution oj(x, y) one chooses and cannot 
be written explicitly in general. As explained above, this phase 
mixing weaken the baroclinic feedback when Pe is large but it 
does not change its intrinsic nature. According to these proper- 
ties, one expect 0^(0) = 1 and (p^x) — > when x — > oo. 

Although cpu can only be determined for a specific vortex 
solution, one can still deduce a few general properties for the 
SBI. For very large thermal diffusivity (very small Pe) one will 
expect the SBI when 



> 1, 



(24) 



(-N 2 )cr 4 

since (p^ ~ 1 in this limit. Although not quantitatively equiva- 
lent, this first limit corresponds to the Pe > 25 threshold of our 
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Fig. 8. Streamline in a vortex undergoing the SB I. A fluid parti- 
cle is accelerated by buoyancy effects on A-B and C-D branches. 
Cooling occurs on B-C and D-A branches (see text). 



simulations. It is moreover similar to the criterion for the onset 
of convection (in the absence of shear) based on the Rayleigh 
number Ra = -N 2 cr 4 /yjd. One finds another instability condi- 
tion due to phase mixing when \i —> 0. Assuming 0^(x) decays 
faster than \jx when x — > oo, one gets a minimum value for fi 
solving 



o(Sct 2 Iia) v 



(-N^cr 4 



which can be calculated once W is explicitly provided. In our 
simulations, this second limit corresponds to the Pe < 16000 
threshold. 

Although not entirely conclusive, this short analytical anal- 
ysis tends to explain why a relatively strong thermal diffusion is 
re quired to get the SBI. This dependancy was also pointed out 
by lPetersen et al.l d2007b) who used both a thermal diffusion and 
a cooling relaxation time in the entropy equation. Such a cool- 
ing time can also be introduced in our analysis in place of the 
thermal diffusion but this does not change the underlying physi- 
cal process. We also note that in this analysis, we have assumed 
that a finite amplitude background vorticity satisfying (TT31) ex- 
isted in the first place, from which we have derived the baro- 
clinic feedback on the same vorticity field. In this sense, this 
analysis describes a n on-linear instability and is d ifferent from 
the linear approach of iJohnson & Gamrniq (l2005al) . In addition, 
we recall that we made the approximation of only including the 
background shear in the advection term on the left hand side of 
equation (IT4b . We remark that the dominant driving term in the 
limit of large [i in (f22l) does not depend on this assumption be- 
cause the heat diffusion term dominates in this limit. 



3.5. Phenomenological picture 

Knowing the basic mechanisms underlying the SBI, one can ten- 
tatively draw a phenomenology of this instability. For this ex- 
ercise, let us consider a single streamline in a vortex subject 
to the SBI (Fig. [8] left). For simplicity, we reduce the trajec- 
tory followed by a fluid particle in the vortex to a rectangle 
(Fig. [8] right). Let us start with a fluid particle on A, moving 
radially inward toward B. This fluid particle is initially in ther- 
mal equilibrium, and we assume the motion from A to B is fast 
enough to be considered approximately adiabatic. As the par- 
ticle moves inward, its temperature and density slowly deviate 
from the background values. If the background entropy profile 
is convectively unstable (condition [6]), the fluid particle moving 
inward is cooler and heavier than the surrounding, and is con- 
sequently subject to an inward acceleration due to gravity (the 
particle "falls"). Once the particle reaches B, it drifts azimuthally 
toward C. On this trajectory, the background density and temper- 
ature is constant. Therefore, the fluid particle gets thermalised 



with the background and, if the cooling is fast enough, it reaches 
C in thermal equilibrium. Between C and D, the same buoyancy 
effect as the one between A and B is observed: as the particle 
moves from C to D, it is always hotter and lighter than the sur- 
rounding, creating a buoyancy force directed radially outward 
on the particle. Finally, a thermalisation episode happens again 
between D and A, closing the loop. 

From this picture, it is evident that fluid particles get accel- 
erated by buoyancy forces on A-B and C-D branches. In the 
end, considering many particles undergoing this buoyancy cycle, 
the vortex structure itself gets amplified, explaining our results. 
Moreover, the role played by cooling or thermal diffusion is ev- 
ident. It the cooling is too fast, fluid particles tend to get ther- 
malised on the A-B and C-D branches, reducing the buoyancy 
forces and the efficiency of the cycle. On the other hand, if no 
cooling is introduced, the particles cannot get thermalised on the 
B-C and D-A branches. In this case, the work due to buoyancy 
forces on the B-A trajectory will be exactly opposite to the work 
done on the C-D trajectory, neutralising the effect on average. 



(25) 4. 2D SBI in compressible flows 



In this situation we consider the 2D subcritical baroclinic insta- 
bility within the framework of a compressible model. We ac- 
cordingy adopt the compressible counterparts of equations (Q])- 
® for an ideal gas with constant ratio of specific heats y in the 
form 



— = --VP - V® - 2ft x u 

Dt p 

-2ClSxe x + V -T v , 
Dc 2 _ (y-\)c 2 Dp 1 
Dt p Dt p 

d t p - -V • pu. 



-i(H + KAc 2 ) 



(26) 
(27) 
(28) 



Here c is the isothermal sound speed which is proportional to 
the square root of the temperature, T v is a viscous stress tensor, 
*K/(7 - 1) is the rate of energy production per unit volume, O is 
a general external potential, and the constant K = \ip where as 
for the incompressible case [i is a thermal diffusivity. 

For simplicity we adopt a model that may be either regarded 
as assuming that there is a prescribed energy production rate *H 
and y is close to unity in (|27t so that the compressional heating 
term oc (y - 1) can be dropped, or replacing the combination 
of the compressional heating term and W by a new specified 
energy production / cooling rate. In either view, as the energy 
production rate is specified, the model is simplified as equation 
(f2Tb becomes an equation for c 2 alone. 

We consider shearing box models of extent L x in the radial 
direction and L y in the azimuthal direction. A characteristic con- 
stant sound speed, Co, defines a natural scale height H = cq/Q,. 
We present simulations using NIRVANA (see section 12.31) for 
a small box with L x = L y /2 = 0.6// and a large box with 
L x = Ly/2 = 1.2//. Thus in terms of the characteristic sound 
speed, when Keplerian shear is the only motion, relative mo- 
tions in the small box are subsonic, whereas supersonic relative 
velocities occur in the big box. We have also considered the two 
resolutions (N x ,N y ) = (144,288) and (N x ,N y ) = (288,576) and 
tests have shown that these give the same results. We have also 
performed simulations with no applied viscosity and with an ap- 
plied viscosity corresponding to a Reynolds number L 2 £llv = 
12500. This was applied as a constant kinematic Navier Stokes 
viscosity but with stress tensor acting on the velocity v being 
the deviation from the background shear rather than u. This is 
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Fig. 9. Time history of the evolution of the enstrophy corre- 
sponding to the perturbations for the small box. The uppermost 
curve is for {ft - 0.25 with no applied viscosity. The next upper- 
most is the corresponding case with an imposed viscosity. The 
lowermost curve is for {ft = 0.025 with applied viscosity and 
the next lowermost curve is for the corresponding case with no 
viscosity 
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Fig. 10. Vorticity map for the small box simulation with applied 
viscosity and {ft = 0.25 at t = 2000. 



not significant for an incompressible simulation but ensures that 
unwanted phenomena such as viscous overstabilty p roduced by 
pertu rbations of the background viscous stress (see iKley et all 
Il993l) are absent. The diffusivity corresponding to the viscosity 
we used is the same magnitude as the magnetic diffusivity used 



Fig. 11. Time history of the evolution of the enstrophy corre- 
sponding to the perturbations for the large box. The uppermost 
curve is for {ft = 0.5 with no applied viscosity. The lower curve 
is for the corresponding case with an imposed viscosity. 



in iFromang et alJ (120071) and as in their case we find that it is 
adequately resolved. 

In order to perform simulations corresponding to the incom- 
pressible ones we need to impose gradients in the state vari- 
ables that produce a non zero Richardson number (appropriate 
for y = 1.) 



i dp a t (P\ 

Ri = — — — In - . 

S 2 pdRdR \p 



(29) 



Because of the shearing box boundary conditions, imposed 
background gradients have to be periodic. Thus we choose an 
initially uniform density p = po and initial profiles for P and 
c 2 of the form P = poc^il - C p sin(2nx/L x )) and c 2 = c^(l - 
C p sm(2nx/L x )), where C p is a constant. For a given value H/L x , 
Ri depends only on C p which was chosen to give a maximum 
value for this quantity of -0.07. Note that in our case, as C p 
is small, we have approximately Ri oc cos 2 (2nx/L x ) in contrast 
to the incompressible simulations where it is constant. Once the 
initial value of fi has been chosen to correspond to a specified 
uniform Peclet number, the external potential O and the energy 
source term *H were chosen such that the choice of state vari- 
ables with Keplerian shear resulted in an equilibrium state and 
then held fixed. For comparison with the incompressible simula- 
tions we adopted Pe = 379 for the simulations reported here. At 
this point we stress that in view of the rather artificial set up, this 
model is clearly far from definitive. However, it can be used to 
show that the phenomena found in the incompressible case are 
also seen in a compressible model but with the addition of the 
effects of density waves in the case of the large box. 

As for the incompressible case, solutions with sustained vor- 
tices were only found for sufficiently large initial velocity pertur- 
bations of the equilibrium state. We adopted an incompressible 
velocity perturbation of the form 

v x = OL y £ll(%7i)){ft sin(27ry /Ly) cos(2nx/L x ) (30) 
and 

= -(3L^/(^7rL x )){ftcos(27ry/L y ) sin(27ix/L x ), (31) 
where {ft is specified amplitude factor. 
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Fig. 12. Vorticity map (left) and surface density map (right) for the large box simulation with applied viscosity and {ft - 0.5 at 
t=190. 
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Fig. 13. Time history of the running means of the the spatially 
averaged value of a for the small box (lower curve) and large 
box (upper curve) simulations with applied viscosity and with 
{ft = 0.5 and {ft = 0.25 respectively. 



4.1. Simulation results 

In Fig. [9] we show the the evolution of the enstrophy associated 
with the perturbations for the small box. Cases with {ft = 0.25 
with and without applied viscosity are shown. These lead to so- 
lutions for which anticyclonic vortices are sustained for long 
times. The case with no viscosity is more active as expected but 
otherwise looks similar to the case with applied viscosity. By 
contrast when the amplitude of the initial perturbation is reduced 
to {ft = 0.025 no sustained vortices are seen. This demonstrates 



that our numerical setup is not s ubject to any linear in stability, 
such as Rossby wave instabilities (lLovelace et al.ll 19991) . The en- 
strophy attains a low level in the case with no applied viscosity 
which is a consequence of a long wavelength linear axisymmet- 
ric disturbance that shows only very weak decay provided by 
numerical viscosity in this run. Fig.[l0]shows a vorticity map for 
the small box simulation with applied viscosity and {ft = 0.25. 
Anticylconic vortices are clearly seen in this case supporting the 
finding from the incompressible runs that a finite amplitude ini- 
tial kick is required to generate them. 

The time history of the evolution of the enstrophy for large 
box simulations with {ft = 0.5 with and without applied vis- 
cosity is shown in Fig. [TT] Corresponding vorticity and surface 
density maps for the case with applied viscosity are shown in 
Fig. [12] Again the inviscid case is more active but nonetheless 
the corresponding maps look very similar. In these cases the an- 
ticyclonic vortices are present as in the small box case but there 
is increased activity from density waves as seen in the surface 
density maps. These waves could be generated by a process 
simi lar to the swing amplifier with vorticity source described 
by iHeinemann & Papaloizoul (2009a b), although the structures 
we observe do not strictly correspond to a small scale turbulent 
flow. The density waves are associated with some outward an- 
gular momentum transport. However, the value of a measured 
from the volume average of the Reynolds stress is always highly 
fluctuating. Accordingly we plot running means as a function of 
time for the small box and large box with applied viscosity in 
Fig. Q3] In the small box case there is a small residual time av- 
erage ~ 10" 4 but in the large box this increases to ~ 3 x 10" 3 . 
This is clearly a consequence of the fact that the small box is 
close to the incompressible regime, whereas the large box being 
effectively larger than a scale height in radial width allows the 
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Fig. 14. Time history of the maxima of the 3 components of the 
velocity field when starting from 2D+3D noise. Once the SBI 
has formed vortices, 3D motions due to parametric instabilities 
appears and balance the 2D instability. 





vortices to become large enough to become significantly more 
effective at exciting density waves. 



5. The SBI in 3D 

The results presented previously were obtained in a 2D setup. 
It is however known that vortices like the one observed in these 
simulations are linearly u nstable to 3D perturba tions due to para- 
metric instabilities (iLesur & Papaloizoull2009l) . The question of 
the survival of these vortices in 3D is therefore of great impor- 
tance, as their destr uction would lead to t he di sappearance of the 
SBI. As shown by ILesur & Papaloizoul (120091) . 3D instabilities 
involve relatively small scales compared to the vortex size when 
the aspect ratioj^ of the vortex is large. This leads to strong nu- 
merical constraints on the resolution one has to use to resolve 
both the 2D SBI and 3D instabilities. To optimise the computa- 
tional power, we have carried out all the 3D simulations using 
our spectral code in the Boussinesq approximation, with a setup 
similar to §[3] and constant stratification. 



5.1. Fiducial simulation 

We first consider a box extended horizontally to mimic a thin 
disc with L = L x = L y = 8 and L z = 1. We set Re = 
1.02 x 10 6 , Pe = 6400 and Ri = -0.01 with a numerical res- 
olution N x x N y x N z = 1024 x 512 x 128 in order to resolve 
the small radial structures due to the 3D instabilities. The hori- 
zontal boundary conditions are identical to the one used in 2D 
and periodic in the vertical direction. Note that we don't include 
any stratification effect in the vertical direction to reduce com- 
putational costs. Initial conditions are to be chosen with care as 
3D random perturbations will generally decay rapidly due to the 
generation of 3D turbulence everywhere in the box. To avoid this 
effect, we choose to start with a large amplitude 2D white noise 
(( Vv^> ~ 0.1) to which we add small 3D perturbations with an 
amplitude set to 1% of the 2D perturbation amplitude. 

We show on FigO the time history of the extrema of the 
velocity field in such a simulation. As expected, vertical motions 
triggered by 3D instabilities appear at t = 200. However, these 



2 The aspect ratio of a vortex is defined by the ratio of the azimuthal 
size to the radial size of the vorticity patch generated by the vortex 



Fig. 16. Time history of v z extrema (top) and of the vortex aspect 
ration (bottom). Once the SBI has formed vortices, 3D motions 
due to parametric instabilities appears and balance the 2D insta- 
bility. 



"secondary" instabilities do not have a destructive effect on the 
SBI. Instead, an almost steady state is reached at t ^ 600. 

Looking at the snapshots in the "saturated state" at t = 800 
(Fig. Q3]) we observe the production of large scale anticy clonic 
vortices, as in the 2D case. However, in the core of these vortices 
(regions of negative vorticity), we also observe the appearance 
of vertical structures. Looking at the vertical component of the 
velocity field, one finds clearly that these vertical structures and 
motions are localised inside the vortex core. Although it is likely 
that the 3D instability observed i n this simulation is re l ated to the 
elliptical instability described by lLesur & Papaloizoul (120091) . we 
can't formally prove this point as the vortic es found in the simu- 
latio n are different from the one studied by lLesur & Papaloizou 
(120091) . 

Despite the presence of 3D turbulence in these vortices, the 
turbulent transport measured in this simulation is directed in- 
ward, with a ^ -3 x 10" 5 . This is to be expected as the 3D insta- 
bilities extract primarily their energy from the vortex structure 
and not from the mean shear. However, allowing for compress- 
ibility will certainly change the turbulent transport as vortices 
will then also produce density waves (see section [4]). 

5.2. Evolution of a single vortex 

To isolate the interplay between the SBI and 3D elliptical insta- 
bilities, we have considered the case of an isolated 2D vortex in 
a 3D setup. We take the same parameters as in the fiducial case, 
but we consider this time a Kida vortex of aspect ratio x - 6 as 
an initial condition. Thanks to these initial conditions, it is possi- 
ble to follow the evolution of this single vortex, and in particular 
measure its aspect ratio as a function of time. This is done with 
a post-processing script, assuming that the vortex core boundary 
is located at cjb = 0.5[mm(co) + max(^)]. We then measure the 
vortex aspect ratio as being the ratio l y /l x of the boundary de- 
fined above. One should note however that this procedure does 
not check that the vortex is still an ellipse. Moreover, it gives 
wrong values for^ if 3D perturbations create strong fluctuations 
of CO. 
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Fig. 15. 3D structure obtained from our fudicial 3D simulation at t = 800. Left: vertical component of the vorticity. Right: vertical 
component of the velocity field. We find that 3D instabilities involving vertical motions and vertical structures are localised inside 
vortex cores. 
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Fig. 17. Vorticity map from a simulation starting with a 2D Kida vortex plus a weak 3D noise. We show a snapshot at t ■ 
t - 750 (middle) and t = 830. After the 3D instability burst, a weak vortex survives and is amplified by the SBI. 
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We show on Fig. [16] the result of such a procedure as well 
as the extrema on the vertical velocity for comparison. Starting 
from x — 6, the first effect of the SBI is to reduce the vortex 
aspect ratio. This is consistent with the Kida vortex solution for 
which 



OJ 

J 



--( — \ 



(32) 



Therefore, the vorticity amplification in the vortex core due to 
the baroclinic feedback leads to a reduction of the aspect ra- 
tio. When x - 4, we note the appearance of strong vertical 
motions. During these "bursts" of 3D turbulence, the aspect ra- 
tio measurement is not reliable anymore. However, the vertical 
motions become ultimately weak, and a weaker vortex appears 
with^f ~ 10. The baroclinic feedback then amplifies the vortex 
and the cycle sta rts again. Interesting l y, the x < 4 Kida vor- 
tices described by lLesur & Papalo izou ( 2009) were shown to be 
strongly unstable due to a "horizontal" instability. The behaviour 
observed in these simulations tends to indicate that the paramet- 
ric modes unstable forx > 4 are not fast enough to take over the 
SBI. However, as x passes through the critical value x - 4, the 
horizontal modes become unstable and strongly damp the vor- 
tex in a few orbits. This interpretation is somewhat confirmed 
by the snapshots of the velocity field (Fig. [17]), where growing 
perturbations are localised inside vortex cores as observed by 
iLesur & Papaloizou! (12009b . 

The "burst" of turbulence observed in the Kida vortex case 
might be a specific property of this vortex solution. Indeed, on 
longer timescales, the spatial distribution of vorticity might be 
modified, leading to a more progressive appearance of 3D insta- 
bilities and a state more similar to the one observed in 15. II could 



be achieved. However, this extreme example clearly demonstrate 
that 3D parametric instabilities do not lead to a total destruction 
of the vortices produced by the SBI. 

6. Conclusions 

In this paper, we have shown that a subcritical baroclinic insta- 
bility was found in shearing boxes. Our simulation s produce vor- 
tices wi th a dynamic similar to the description of iPetersen et al.l 
(l2007al) and lPetersen et al.1 (l2QQ7bh . We find that the conditions 
required for the SBI are (1) a negative Richardson number (or 
equivalently a disc unstable for the radial Schwarzschild crite- 
rion) (2) a non negligible thermal diffusion or a cooling function 
and (3) a finite amplitude initial perturbation, as the instability is 
subcriticaQ. 

When computed in compressible shearing boxes, the vortices 
sustained by the SBI produce density waves which could lead 
to a weak outward a ngular momentum transport {a ~ 10~ 3 ), a 
process suggested by lJohnson & Gammid d200 5b). However, we 
would like to stress that this transport cannot be approximated 
by a turbulent viscosity, as the transport due to these waves is 
certainly a non local process. Several uncertainties remain in the 
compressible stratified shearing box model. In particular, the im- 
posed temperature profile that is subsequently maintained by a 
heating source is somewhat artificial. This occurs while angu- 
lar momentum transport causes a significant density redistribu- 
tion causing the local model to deviate from the original form. 
Clearly one should therefore consider global compressible sim- 
ulations with a realistic thermodynamic treatment to draw any 

3 Note tha t point (1) and (2) were already mentioned by 
Petersen et al. ( 2007b), although in a somewhat different form. 
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firm conclusion on the long term evolution and the angular mo- 
mentum transport aspect of the SBI. 

Although the vortices produced by the SBI are anticyclones, 
the precise pressure structure insi de these vortices is not as sim- 
ple as the description given by iBarge & Sommerial (Il995l) . In 
particular, since the vorticity is of the order of the local ro- 
tation rate, vortices are not necessarily in geostrophic equilib- 
rium and pressure minima can be found in the centre of strong 
enough vortices. Moreover, these vortices interact with each 
other, which leads to complex streamline structures and vortex 
merging episodes. This leads us to co nclude that the particle con- 
centration mechanism proposed by IBarge & Sommerial d 19951) 
might not work in its present form for SBI vortices. 

In 3D, vortices are found to be unstable to parametric insta- 
bilities. These instabilities generates random gas motions ("tur- 
bulence") in vortex cores, but they generally don't lead to a 
proper destruction of the vortex structure itself. More impor- 
tantly, the presence of a weak hydrodynamic turbulence in vortex 
cores will lead to a diffusion of dust and boulders, balancing the 
concentration effect described above. In the end, dust concentra- 
tion inside these vortices might not be large enough to trigger a 
gravitational instability and collapse to form planetesimals. 

Given the instabi lity criterion detailed above , one can ten- 
tatively explain why iJohnson & Gammiel (120061) failed to find 
the SBI in their simulation. First, we would like to stress 
that our compressible simulations are not more resolved than 
theirs. Since these simulations use similar numerical schemes, 
the argument of a to o small Reynolds number proposed by 
Petersen ^et al.l (l200"7b|) see ms dubious. We note however that 
Johnson & Gammie ( 2006) did not include two other important 
ingredients: a thermal diffusivity and finite amplitude initial per- 
turbations. As shown in this paper, if one of these ingredient is 
missing, the flow never exhibits the SBI and it gets back to a 
laminar state rapidly. Note that even if thi s argu ment explains 
the negative result of Johnson & Gammie ( 2006|), it also indi- 
cates that the SBI should be abs e nt fro m the simulations pre- 
sented by Kla hr & Bodenheimerl (120031) . as no thermal diffu- 
sion nor cooling function was used (the initial conditions were 
not explicitly mentioned). This s uggests that either the globa l 
baroclinic instability described by Klahr & Bodenheimerl (120031) 
and the SBI are two separate instabilities, or strong numerical 
artefacts have produced an artificially large thermal diffusion in 
Klahr & Bodenheimer ( 2003) simulations. 

In the limit of a very small kinematic viscosity, the scale at 
which the instability appears has to be related to the diffusion 
length scale (fi/S) l/2 . In a disc, one would therefore expect the 
instability around that scale, provided that the entropy profile 
has the right slope (condition 1). However, as the instability pro- 
duces enstrophy, vortices are expected to grow, until they reach a 
size of the order of the disc thickness where compressible effects 
balance the SBI source term. We conclude from this argument 
that even if the SBI itself is found at small scale (e.g. because of 
a small thermal diffusivity), the end products of this instability 
are large scale vortices, with a size of the order of a few disc 
scale heights. 

We would like to stress that the present work constitutes 
only a proof of concept for the subcritical baroclinic instabil- 
ity. To claim that this instability actually exists in discs, one has 
to study the disc thermodynamic in a self consistent manner, a 
task which is well beyond the scope of this paper. Moreover, it 
is not possible at this stage to state that the SBI is a solution 
to the angular momentum transport problem in discs. Although 
this instability creates some transport through the generation of 
waves, this transport is weak and large uncertainties remain on 



its exact value. Finally, t he coexistence of the SBI and the MRI 
dBalbus & Hawleyl lT991) is for the moment unclear. In the pres- 
ence of magnetic fields, vortices produced by the SBI might 
beco me strongly unstable be cause of magnetoelliptic instabili- 
ties (Mi zerski & Bajerl [20091) . Although the nonlinear outcome 
of these instabilities is unknown, one can already suspect that 
vortices produced by the SBI will be weaken in the presence of 
magnetic fields. 
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